Universal microstructure and mechanical stability of jammed packings 
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Jammed packings' mechanical properties depend sensitively on their detailed local structure. 
Here we provide a complete characterization of the pair correlation close to contact and of the force 
distribution of jammed frictionless spheres. In particular we discover a set of new scaling relations 
that connect the behavior of particles bearing small forces and those bearing no force but that are 
almost in contact. By performing systematic investigations for spatial dimensions d=3-10, in a 
wide density range and using different preparation protocols, we show that these scalings are indeed 
universal. We therefore establish clear milestones for the emergence of a complete microscopic 
theory of jamming. This description is also crucial for high-precision force experiments in granular 
systems. 



Introduction - The jamming phenomenon is ubiqui- 
tous - candies [1], coal [2], and colloids [3] all can jam, 
but its microscopic universality remains debated, even 
for the most ideal of systems. Like any other phase tran- 
sition, the jamming transition can be approached from 
the unjammed phase, e.g. by compressing hard spheres 
(HS) [4], or from the jammed phase, e.g. by minimizing 
the energy of soft spheres (SS) [5]. Yet these two com- 
plementary approaches have mostly been developed in- 
dependently from each other (see [6] for HS and [7, 8] for 
SS). Unlike standard phase transitions, however, the jam- 
ming transition is a non-equilibrium phenomenon that 
happens deep inside the glass phase [9, 10], and there- 
fore different protocols generate different packings, which 
may result to conflicting observations. Indeed, all agree 
that marginally stable packings of frictionless spheres 
average 2d force-bearing contacts per particle [8], but 
jammed packings' density [6, 11-14], parts of their mi- 
crostructure [6, 15, 16], as well as their given name [6, 17] 
are contentious. Although the jamming "j" point was 
proposed to be unique in the thermodynamic limit [5, 18], 
there is a growing consensus that jamming occurs over 
a range of "j" points [6, 7, 9, 13, 14, 17]. Yet vari- 
ous physical origins have been attributed to the jamming 
density variation, including structural correlations in the 
initial configuration [7], and the presence of small crys- 
talline regions only detectable by subtle order metrics 
whose minimization should result in a single "maximally 
random jammed" state [6, 17]. Others have proposed 
the intrinsic existence of a range of densities over which 
packings with an identically disordered structure could 
be found [9, 13, 14]. A power-law growth of the num- 
ber of almost-touching particles near jamming has also 
been identified numerically, but different exponents have 
been found for HS [4] and SS [15]. If there is microscopic 
universality, it has yet to fully emerge. 

In this letter we bring a different point of view to the 
problem by systematically investigating how the jam- 
ming limit is approached from both sides of the transition 



and by varying the dimensionality of space from d=3 to 
10. This approach allows us to obtain a series of impor- 
tant results. (1) Increasing d > 4 suppresses crystalliza- 
tion [19, 20] and the "spurious" contribution of "rattlers" . 
We can thus show that random jammed packings of 
monodisperse spheres with identical near-contact struc- 
tural properties can be obtained over a range of densities 
(thus confirming results in d=3, 4 [4, 14, 17, 19, 21]), 
and that this range broadens with increasing d. (2) We 
confirm an earlier suggestion that two exponents a and 
0, corresponding to different physical regimes, control 
the mechanical stability of jammed packings [22]. The 
first describes the "quasi-contact" regime in which par- 
ticles are separated by very small gaps h, whose number 
scales as h~ a for small h; the second describes the tail 
of the "contact" regime, where the number of particles 
bearing a small force / scales as f e . (3) We also pro- 
vide a complete characterization of the microstructure of 
jammed packings. We show that matching the two above 
regimes provides scaling relations between the exponents 
and non-trivial scaling functions. We thus conclude that 
the mechanical stability of jammed packings is related to 
their very complex contact microstructure. (4) We find 
these results to be universal in the sense that they are ro- 
bust to changes in preparation protocol, packing density, 
and, in particular, spatial dimension. 

The observation that jammed packings' properties are 
independent of d suggests that a mean-field theory should 
be able to capture the jamming phenomenology [23, 24]. 
One such treatment, the Gaussian replica theory (G- 
RT) [10, 13], unifies the description of the glass transition 
and of jamming by exploiting an analogy with discrete 
random optimization problems [9, 25]. In this treatment, 
the HS and SS approaches to jamming are unified un- 
der the assumption that jammed states are the infinite 
pressure (for HS) or zero temperature (for SS) limit of 
long-lived metastable glassy states [10, 13]. The theory 
predicts a growing jamming density range with d [13], 
the existence of scaling relations for energy and pressure 
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relating the two sides of the jamming transition [10], and 
makes structural scaling predictions that are remarkably 
satisfied at short distances [10, 13]. Yet we show here 
that (5) G-RT completely fails to describe the structural 
regime that controls jammed packings' mechanical sta- 
bility. Our results (l)-(5) will thus guide both theory 
and experiments (through high-precision force measure- 
ments [26]) towards a better understanding of the jam- 
ming transition. 

Packing Generation - We consider a system of N > 
8000 identical spherical particles of diameter a in a fixed 
volume V , under periodic boundary conditions. The 
packing fraction tp = NVd(o~/2)/V, where Vd(r) is the 
volume of a d-dimensional sphere of radius r, measures 
the fraction of space occupied by particles. Jammed 
packings are prepared using two different numerical pro- 
tocols (see Appendix for details and reduced units defini- 
tions), (i) Approaching jamming from densities below it 
by Lubachcvsky-Stillingcr (LS) compressions of HS un- 
dergoing Newtonian dynamics while a grows at a fixed 
rate j=& [4]. The compression, which is tuned to prevent 
crystallization [20, 27], stops when particles are very near 
contact, defining the packing fraction tp^ at which the 
HS reduced pressure becomes infinite, (ii) Approaching 
jamming from densities above or below it by minimiz- 
ing the energy E of a random configuration of harmonic 
SS. Initial bounds er_ and er + that bracket jamming are 
evolved itcratively by choosing an intermediate value cr m 
and minimizing the energy of the current configuration 
at cr + (procedure from above) or at ct_ (procedure from 
below). The final jammed configurations at the onset of 
E 7^ have tp\. from above and tp\ from below. From 
above, the energy vanishes with e = E/N ~ Atp 2 and 
the static pressure P ~ Aip, where Aip is the distance 
from jamming [5]. 

We find the initial a± to have no measurable effect 
on tp\. We formally define p™ m = min CT± tp\{a±), but 
any reasonable a± results in the same final density. By 
contrast, tp* strongly depends on a + (Fig. 1), but is 
also independent of er_. We therefore define (^™ ax = 
max ff± tp*(a±). A practical way of constructing both 
cp™ ln and <^>™ ax is to run the energy minimization (respec- 
tively from below and from above) starting from er_ = 
and cr + large enough to saturate tp* to its maximum. 
Intermediate packing fractions can then be obtained by 
reducing a + (Fig. 1). By varying a± in protocol (ii) 
we can thus construct packings over a density interval 
j-^mm^maxj fo&t roughly corresponds in protocol (i) to 
[ipp~ , tpj + ] with 7_ ps 3x 10~ 2 and 7+ ps 3x 10~ 4 (larger 7 
generate mechanically unstable packings). The resulting 
density range is remarkably found to grow steadily from 
about 2% in d = 3 to nearly 10% in d = 11 (Fig. 1). We 
therefore confirm the similar observation made for d = 3 
binary mixtures [14], where the limited available density 
range and the subtle crystal order had left some room for 
debate [6]. Note that this range is achieved by only im- 
plementing procedures that compact liquid congurations. 
Ref. [21] has shown that enlarging the space of procedures 
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FIG. 1. The extrapolated jamming density tpp following 
the protocol described in Ref. [27] is extended to higher d 
(solid line and crosses), and compared with the G-RT pre- 
diction for tpGCP (dashed line), (top inset) The range of 
jamming densities tp], (squares) is compared to <p™ a * (cir- 
cles) and ipf ln (triangles). Note that tpf^ ~ ^ =3xl0_4 and 

<P? in ~ tPp =ixl0 ~ 2 ■ (bottom inset) The d=3 increase of ipi 
with a+, in terms of the initial effective packing fraction. 

enlarges the range of jammed packings, but the resulting 
packings likely have a different microstructure. 

The similarity between the jamming density results of 
the two protocols suggests an underlying physical connec- 
tion between them. G-RT indeed predicts that packings 
exist over a finite packing fraction range, whose upper 
limit is the "glass close packing" </?gcp [13]. By analogy 
with random combinatorial optimization problems [25], 
the densest packing at pgcp is conjectured to require a 
time ~ exp(N a ) to generate, the exponent being possi- 
bly a ps (0? — l)/d, based on a nucleation analysis. The 
maximal density that can be reached by the protocols 
above, which both run in polynomial time in N, should 
therefore be strictly smaller than </?gcp- Figure 1 shows 
it to be the case for all d, in agreement with G-RT. 

Scaling functions - To determine the universal struc- 
ture of disordered jammed structures, we consider the 
pair correlation function g(r) — (pN)^ 1 <H r + r i ~ 

r_j ) ) , which is the only relevant structural correlation in 
high d fluids [28]. For numerical convenience, we com- 
pute the cumulative structure function 

Z(r) =pS d _i I dss d - l g{ S ) , (1) 
Jo 

where S d —i is the surface of a d dimensional sphere of 
unit radius. The function Z[r) thus provides the average 
number of neighbors within r of a given particle. Rattlers 
are first excluded from the analysis (Appendix), but we 
come back to this point below. 

For both protocols, Z(r) jumps from to a plateau 
at z on a scale proportional to the distance to jamming 



3 



Z(r 




FIG. 2. Schematic of Z(r) when approaching jamming from 
above (a) with protocol (ii), and from below (b) with proto- 
col (i). Three distinct scaling regimes can be identified. The 
first regime is related to the growth of Z(r) from to z. It 
corresponds to interparticle gaps h = \r — a\ ~ Atp, and hence 
to particles that are in contact when Atp — > 0. The last regime 
corresponds to gaps h that remain finite for Atp — > 0. These 
particles remain separated at jamming, but these small gaps, 
Z(r) — zee h 1 ^" , form a "quasi-contact" regime. The inter- 
mediate regime corresponds to gaps h ~ Atp 11 . It matches the 
two other regimes and disappears when Atp — > 0. 



Atp, where z is the isostatic average number of contacts 
2d plus a correction z — 2d oc Atp** (Fig. 2). For HS, 
we find C=0.36(2), while C=0.53(3) for SS (Fig. 3) [5]. 
The approach to the isostatic plateau is characterized by 
a long power-law tail with exponents 0=0.28(3) for HS 
and 0=0.42(2) for SS, but the exponent is independent 
of d for a given model. The plateau is extended by a 
second power-law regime that corresponds to particles in 
"quasi-contact" , carrying no force at jamming. We find 
that in this regime the scaling is the same for both pro- 
tocols, growing as Z{r) — z oc (r — o) l ~ a with a universal 
exponent a = 0.42(2) until it reaches the trivial large r 
regime. Interestingly, the two power-law regimes can be 
matched by a scaling function T~i± 1 which defines an ad- 
ditional intermediate regime. This intermediate regime 
shrinks to a point at jamming, but smoothly crosses over 
from one power-law regime to the other at finite Atp. 
Consistency therefore sets clear scaling requirements for 
the different regimes (see Appendix for scaling analysis) 
as detailed in Fig. 2, and verified in Fig. 3. 

Force distribution and mechanical stability - The 
consequences of these universal scaling relations on me- 
chanical properties can be gleaned from the probability 



distribution of inter-particle forces /. Here again, we 
consider the cumulative distribution G(f) = J P{f')df 
rather than the pair force distribution P(f), for numeri- 
cal convenience. 

For HS approaching jamming, the average force / oc p. 
In the contact regime the force and distance distribu- 
tions are also related through a Laplace transform (Ap- 
pendix) [4]. The low- force distribution is thus consistent 
with G{f) oc f 1+e and 9 = 0.28(3). For SS approaching 
jamming from above, the pair potential sets the relation 
between the force and the pair distributions [29] (Ap- 
pendix). Here again, the low- force tail is consistent with 
9 = 0.42(2). For both protocols, however, the regime in- 
termediate between contacts and quasi-contacts results 
in deviations from this power-law decay at very weak 
forces away from jamming. 

The large force regime has been thoroughly stud- 
ied [4, 5, 18, 29-32], but the weak force distribution is 
much less well characterized. Yet it has been proposed by 
Wyart [22] that a > 1/(2 + 9) is required for mechanical 
stability. Both the SS values (a = 0.39(1), 9 = 0.42(2)) 
and the HS ones (a = 0.42(2), 9 = 0.28(3)), however, in- 
dicate a slight violation of this condition. A generalized 
stability condition of the form a > (1 — <5/2)/(2 + 6* — 
5/2) [22] is consistent with our findings for S > 0.2, but a 
direct test of this extended relation is beyond the scope 
of the current analysis. 

Rattlers - Rattlers, i.e., particles with no mechani- 
cal contacts, must be considered before concluding that 
the dimensional and protocol robustness of these results 
strongly support a universal microscopic description of 
jamming. Because their fraction rapidly decreases with 
increasing d (Appendix) [19], and their structural contri- 
bution is clearly distinct from that of the other particles 
when Atp — > 0, it is reasonable to remove them from 
the analysis. Rattlers indeed play essentially no role in 
the scaling regimes in high d, while in low d, their inclu- 
sion introduces noise in Z(r) and G(f) that obscures the 
scaling relations, which may explain why a ss 0.5 was 
obtained in Ref. [15]. Removing the rattlers reveals the 
robust relationship between microstructure and mechan- 
ical properties, in support of jamming having a critical 
dimension d c = 2 [24]. 

Comparison with microscopic theory - G-RT, the 
only available first-principle theory of jammed packings, 
provides predictions for the contact regime scaling func- 
tion Z±(x) [10, 13] (Appendix). We find the form of 
Z±(x) to be extremely accurate when x is of order 1, 
but G-RT fails to capture the ensuing power-law regimes 
(Fig. 3). G-RT indeed predicts an exponent 9 = for 
both protocols, and completely misses the power-law di- 
vergence related to a, predicting a = 0. A similar devi- 
ation is observed at weak forces. We attribute these dis- 
crepancies to the Gaussian assumption for the cage form 
of G-RT, which has recently been found to be erroneous 
in dense disordered fluids [33, 34]. This non-Gaussian 
structure also naturally suggests a microscopic explana- 
tion for the breakdown of the normal-mode decomposi- 
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FIG. 3. Scaling of the cumulative structure function Z(r) and the cumulative force distribution G(f) in d = 3 upon approaching 
jamming from above (a) by SS energy minimization (where e oc |A</?| 2 — s> and e = \J ed/2), and from below (b) by HS 
compression (where p oc |A(/3| _1 — s> oo). (la) For diminishing e, the height of the plateau (inset) converges to the isostatic 
value with £=0.53(3). (2a) The small r < a regime shows the "contact" scaling function Z+ (x), which agrees well with the 
G-RT prediction (red line). (3a) Rescaling Z(r) using fj, = (1 + S)/(2 + 9 — a) and v = a/i highlights the behavior of the 
scaling function \H-(x) — "H_(l)| ~ 1.2s (brown line) along with the 9 = 0.42(2) (red line) and the a = 0.39(1) (blue line) 
power-law regimes. (4a) G(f), with power-law tail exponent 9 = 0.42(2) (dashed line), (lb) For increasing p, Z(r) grows 
on an earlier scale r — a ~ p^ 1 to a plateau at the isostatic value, whose height (inset) decays with £ = 0.36(1). (2b) The 
small r — a regime shows the "contacts" scaling function Z-(x), which agrees well with the G-RT prediction (red line). (3b) 
Rescaling Z(r) using fj, = (l + O)/(2 + — a) and v — ctfi highlights the behavior of the scaling function \"H- (x) — T-L-{1) \ ~ &x 
(brown line) along with the 6 = 0.28(3) (red line) and a = 0.42(2) (blue line) power-law regimes. (4b) G(f), with power-law 
tail exponent 9 — 0.28(3) (dashed line), compared with the G-RT prediction (solid line). 



tion of jammed states [35, 36]. Including a non-Gaussian 
cage to RT ought to provide a better mean-field under- 
standing of the jamming phenomenology. 

Conclusions - Our results show that the jamming 
terminology controversy should be resolved by replac- 
ing the j-point [5] with the j-line [9, 13], and by distin- 
guishing a range of maximally random jammed packings 
from their partially crystallized counterparts [6, 17, 21]. 
They also reveal that the contacts' complex microstruc- 
ture in jammed packings is characterized by universal, 
well-defined scaling regimes and by their corresponding 
scaling functions. We give precise numerical predictions 
for the scaling exponents, and show that the scaling func- 
tions are related to the force probability distribution. 
These specific predictions can be tested in soft matter 
and granular experiments. A preliminary investigation 
indeed examined the scaling of the peak of the pair cor- 
relation function [3], but our comprehensive predictions 



can help experimentalists access the full scaling of Z(r) 
and G(f). This feat should be possible once a force res- 
olution of ~ 5%f is experimentally achieved [26]. 

Finally, it is worth noting that the present study was 
limited to zero temperature T in the sense that no ther- 
mal motion is allowed in SS and that for HS the energy 
interaction scale is infinite compared to T. At finite T, 
the jamming transition is blurred [10], but vestiges of 
the scaling relations should remain visible [3]. Future 
work will detail how temperature and its associated an- 
harmonicitics affect the T = scaling relations identified 
here [10, 35]. 
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Appendix A: Numerical simulations 

Molecular dynamics simulations of HS and SS energy 
minimizations of iV=2 18 spheres in d=3, iV=8000 in 
4 < d < 9, and N=2 x 2 14 in d=10 are performed un- 
der periodic boundary conditions. For d = 3, the choice 
of a very large system is motivated by the need of re- 
ducing the statistical noise in the intermediate scaling 
regime for Fig. 3 of the main text. The d > 3 system 
sizes chosen ensure that even when the system is at its 
densest the box edge remains larger than 2a, which pre- 
vents a particle from ever having two direct contacts with 
another one. There are strong reasons to believe that al- 
though relatively small these N nonetheless provide a 
reliable approximation of bulk behavior. First, with in- 
creasing d the largest diagonals of the simulation box are 
\fd larger than the box edge. Second, correlations of the 
fluid structure are expected to decrease very quickly with 
increasing d [19, 28], and correspondingly finite-size ef- 
fects are reduced. The validity of these rationalizations, 
which are consistent with the decorrelation property of 
high d sphere packings [6] , have been satisfactorily tested 
in d = 8 in Ref. [27]. 



1. Hard Sphere LS Compressions 

Event-driven HS simulations are performed at ther- 
mostated inverse temperature j3 for spheres of mass m, 
starting from random configurations in the limit a — > 0. 
Time has units of y/ (3ma 2 [19, 27], but all dimensional 
quantities arc expressed in units such that /3 = 1 and 
m = 1. The reduced pressure p = PP/p, with p = N/V 
for pressure P, diverges at the jamming packing fraction 
ipl as p ~ |A(p| _1 with Aip = ip — (pj [4]. Crystalliza- 
tion in d > 3 is strongly suppressed, so access to deeply 
supersaturated starting configurations can be attained 
via the slow growth rate 7 = a = 3 x 10 4 [19, 20]. In 
d = 3, where the crystallization of monodisperse hard 
spheres is relatively rapid for moderately high packing 
fractions, 7 = 10 -2 is employed up to p = 10 3 , but the 
slow compression rate is used afterwards. The force be- 
tween particles in high p configurations is measured from 
the rate of momentum exchange between pairs of parti- 
cles in simulations with 7 = 0. These measurements are 
made over at least 10 4 collisions per particle. 



2. Soft Sphere Energy Minimizations 

The harmonic SS energy is E(X, a) = Yli>j v (\^i ~ 0*1) 
with X = {ri} and v(r) = e(a — r) 2 9(a — r). Units are 
chosen such that e = 1. We start from random sphere 
configurations X + = A_ and use cr_ and <r + that bracket 
the jamming point, i.e., -E(A_,cr_) = 0, E(X + , a + ) > 0. 
Jamming is identified as the onset of non-zero energy, it- 
crativcly determined using a bisection method. At each 



iteration, an intermediate value a m is chosen, and the 
energy at cr m is minimized via conjugate-gradient (CG) 
minimization starting from either X + (from above) or 
A_ (from below). The configuration obtained after min- 
imization X m then substitutes (A + ,er + ) = (X m ,a m ) if 
E(X m , a m ) > 0, or <r_) = (X m , <7 m ) i£E(X m , a m ) = 
0. The procedure stops when either of two conditions is 
satisfied: the energy per particle E(X + ,a + )/N falls be- 
low 10 -20 , corresponding to typical overlaps of the order 
of 10~ 10 ; or the change in E(X + ,(T + ) from one mini- 
mization step to the next is less than the bound set by 
double precision arithmetic, i.e., lCP 8 ^. The final value 
of ct_ defines the final packing fraction p\ in the proce- 
dure from below and a + that of in the procedure from 
above. 



3. Rattler Analysis 

The rattlers are self-consistently determined by iden- 
tifying the number of particles with fewer than d + 1 
neighbors within a distance cutoff for the smallest \A<p\ 
obtained by a given approach to jamming. In HS, a cut- 
off of tr(l + 100/p) is used, which roughly corresponds to 
a cutoff of a(l — e/100) for SS, where the scaling vari- 
able e = yj e d/2. This cutoff slightly overestimates the 
number of rattlers at jamming, but the rapid diminution 
of the fraction of rattlers with d (Supplementary Fig. 7) 
guarantees the robustness of the results. 

4. Extraction of the jamming density 

Following [27], for each compression run at fixed 7 we 
perform a linear fit of the line 1/p vs ip with p > 50. 
The point where the linear fit vanishes, indicating infinite 
pressure, defines tp^ for this given run. Next, for a fixed 
dimension, we fit (pj = tpj~^° + Ayf^j to extrapolate the 
jamming density at 7 — > 0. 

For the energy minimization protocol, when approach- 
ing jamming from above, the jamming density can be 
obtained by accurately fitting the energy data with e = 
eo[(<p — vt)/'/ 5 t] 2 (inset in Supplementary Fig. 4). We 
focus in particular on minimization runs performed with 
initial <r_ = and er+ — > 00, which in practice correspond 
to taking the largest <t+ at which no variation of is 
detected. The associated prefactor e can then be com- 
pared with the prediction from G-RT. A good agreement 
is obtained when the correction discussed in [10, Section 
VLB] is taken into account (Supplementary Fig. 4). 

Note that pj~^° and (p™ ax arc quite close to each 
other. They indicate the best packing density that can 
be reached using our two different compression protocols. 
Note that according to G-RT both should be smaller than 
the maximal packing density of glassy states, called glass 
close packing <p>gcp- According to the theory, it is very 
unlikely that packings at v?gcp can be produced in poly- 
nomial time, hence it is expected that both ^ _>0 and 
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Supplementary Figure 4. The prefactor of energy e = 
eo[(y> — >fii)/<pi] 2 from numerics (black curve) and from G- 
RT (red curve, bare data; green curve, corrected data, see 
Sec. C2 and [10, Section VLB]), (inset) The scaling e — 
eo[(<£ — ipi)/(pi] 2 upon approaching <p% from above for d — 3, 
with cr_ = and o+ — > oo. Points are numerical data and 
the line is a quadratic fit used to extract eo (the value of (pi 
is obtained by imposing visually the best alignment of the 
data). 
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Appendix B: Scaling functions 

1. Structure Scaling Analysis 

When approaching jamming with protocol (i), Z(r) = 
for r < er, and for r > a, p parametrizes the scaling 
function for Z(r > cr). A first scaling regime i — a ~ 
sees Z(r) grow from to the average number of "con- 
tacts" z as 

Z(r) = zZ- [(r - a)p/a] (Bl) 

with Z-{x) ~ 1 — Cx" 1 " 9 when x — > oo for a constant 
C [4]. Force-bearing contacts are only observed at jam- 
ming proper, but their signature develops asymptotically. 
A second regime for finite r — a has 

Z(r) = z + C'{r-a) 1 ~ a , (B2) 

where C is a constant. At jamming, these nearly touch- 
ing "quasi-contacts" carry no force. For large r, a trivial 
regime develops independently of | A(p\ (Sec. B 3). Match- 
ing the first two scaling regimes implies the existence of 
an additional intermediate regime "H_ for r — a ~ p~ t* 

Z{r) = z + p v -»U- [(r - tr)fP /a] (B3) 

with fi < 1 and v < fi. Consistency then requires that 
T~L-(x — > 0) oc — x^ 1 ^ 9 and H-(x — > oo) oc x 1_Q with 
scaling relations v = a/j, and /x = (1 + 6)/ (2 + 6 — a). 
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Supplementary Figure 5. Cumulative force distribution G(f) 
in d=3-10 (a) for HS with 9 = 0.28(3) and (b) for SS with 
6 — 0.42(2). The force distribution in higher d is essentially 
the same as in d = 3 and the high force behavior agrees 
equivalently well with the G-RT predictions. The exponents 
extracted from the small force tail are also numerically in- 
distinguishable, (inset) Test of Eq. (B7) from the numerical 
results for Z- (x) (points) and plugging the numerical G(f) 
in Eq. (B7) (solid line) 



When approaching jamming with protocol (ii) from 
above, the remaining overlaps provide a scaling variable 
e = \Jed/2 oc \Aip\ [10]. In spite of the very differ- 
ent preparation protocol, similar structural regimes are 
identified. In the r < a contact regime, Z(r) grows from 
to z by a universal scaling function 

Z(r) = zZ+[(a - r)e- x /a] (B4) 

with Z + (x) ~ 1 — C"x 1+e when x — >• for a constant 
C" . For r > a, here again 

Z(r) = z + C (r - a) 1 -" , (B5) 

hence the two regimes must be matched by an interme- 
diate scaling function "H+ for r — a ~ e M 

Z(r) = z + e^ v U+ [(r - ^e^/a] (B6) 

with \i > 1 and v < 1. Consistency here requires that 
H+(x — > — oo) oc — |x| 1+e and 'H+(x — > oo) oc x x ~ a with 
v = a/i and fj, = (1 + 9)/(a + 6). 

In both cases, the exponents a and 9 are determined by 
collapsing the numerical results using this scaling form, 
which was repeated for different systems. This constraint 
leaves a relatively small uncertainty on the final value, 
which provides the error bar. 

2. Force Scaling Analysis 

When HS approach jamming, the cumulative force dis- 
tribution G(f) approaches a scaling function defined by 
G{yf) — > Q~{y)- The relation between the scaling func- 
tions Z_ (x) and G- (y) suggested in Ref . [4] 

I* oo 

Z_(x) = l-x dy (y) e~ xv (B7) 
Jo 
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is verified in Supplementary Fig. 5. It follows that if 
Z_{x) ~ 1 - Cx- 1 - 9 for x -> oo, then Q_(y) ~ y 1+e for 
2/ — > 0, which is here observed for all d > 3 (Supplemen- 
tary Fig. 5). When SS approach jamming from above, 
the interaction potential gives / = 2(er — r) for < r < a 
and zero otherwise, so G(j) = 1 — ^-^nj^ ■ In the jam- 
ming limit G(2yea) —> G+{y) = l—Z + (y), and therefore 
G+(y) ~ y 1+6 \ as in the previous case. This behavior is 
here observed in all d > 3 (Supplementary Fig. 5). 



d = 3, and because different observables are considered, 
additional results are here presented. They are reported 
in this section and incidentally provide a somehow sim- 
plified derivation of the results of Ref. [10]. Nonethe- 
less, reading this section requires a detailed knowledge 
of Refs [10, 13], so the reader who is not interested in 
the theoretical details can safely skip it. Note that as in 
the main text, this section uses reduced units e = 1 and 
a = 1. 



3. High d Structure 

In the contact regime r — a ~ A(p, G-RT predicts 
that scaling functions Z± T should describe the growth of 
Z(r) from to the isostatic value, as given in Eqs. (CIO) 
and (CI 7). Both results are tested in Supplementary 
Fig. 6 for c?=3-10. The collapse is remarkably good for 
all x. The agreement with the G-RT scaling form is 
also remarkable for small x, but start to deviate from 
the theoretical prediction when Z(r) approaches the iso- 
static z rj 2d plateau. The type of deviation is different 
from each protocol, but is similar from one dimension 
to the next for a given protocol. For larger x in the 
near-contact region, the two protocols robustly produce 
the same power-law growth, Z(r) ^ (r — ct) 1_q with 
a ~ 0.40(1) (Supplementary Fig. 6), which is not pre- 
dicted by G-RT. Because the number of rattlers vanishes 
with dimension neither of these phenomena can be as- 
cribed to their presence. But because G-RT predictions 
rely on the individual cages to be Gaussian, which pre- 
sumably they are not [14, 33, 34], it is natural to ascribe 
the discrepancy to the breakdown of that assumption. 

At very large distances, the pair correlation function of 
any disordered systems trivially has g(r ^> a) = 1, which 
corresponds to Z(r 3> 1) ~ 2 d ip[(r/a) d — 1]. Unsurpris- 
ingly this scaling form captures well the behavior of Z(r) 
for both protocols and all d at large r, but the range of 
validity also extends with d (Supplementary Fig. 6). In 
order to quantify this effect, we fit the curves of Z(r) for 
r > a using the form 

Z(r) = C'(r - a) 1 '" + 2 d ip[(r/a) d - 1)] . (B8) 

When d grows, the region where the second term is much 
bigger than the first is 

C'(r - a) 1 -" < 2 d tp{{r/a) d - 1) ~ d2 d tp(r - a)/a (B9) 

hence (r-a) > [a C / (d2 d tp)] 1 / a . The fitted values of C 
indicate that the crossover point indeed decreases slowly 
with d. 

Appendix C: Replica theory calculations 

The predictions of G-RT presented in this work are 
based on earlier results [10, 13]. Yet because the calcula- 
tions in Ref. [10] have only been explicitly carried out for 



1. General expressions 

The approximation scheme used is based on [10, 
Eq. (22) and (23)], which give the replicated free entropy 
separated between the harmonic and the liquid contribu- 
tions 

S(m, A; T, tp) = S h (m, A) + S Vlq (T/m, tp) 

+ 2 d - 1 V y^&)G(m,A;T) ,with 

pOO 

G(m, A;T) = d dr r^ 1 [q(A, T; r) m - e - pmvi - r) ] 
Jq 

(CI) 

for m replicas at temperature T, in a Gaussian 
cage of variance 2A. The function q(A, T; r) = 

J d d r' l2 A{f')e~ fSv( ^ f: " {) is defined in [10, just after 
Eq. (16)] where 72^4 is a normalized and centered Gaus- 
sian of variance 2A, and y^ is the HS cavity function. 
Introducing bipolar coordinates, as in [13, Appendix 
C.2.a], we obtain the generalization of [13, Eq. (C16)] 
to the soft sphere potential v(r) = (1 — r) 2 8(l — r) 

q (A,T;r)= Tdue-^U^^^, 

Jo \r/ V47rA ( C2 ) 

_ n / ru f ru\ 
X e TK \ ir—r Id-2 [-—) . 

V A — \2A) 

The above equations (CI) and (C2) are the starting 
point of all the needed replica calculations for our anal- 
ysis, and because we focus on the "jamming limit" of 
these equations, we take T — > with r = T/m and 
a = A/m held constant [10]. In Ref. [10] this limit was 
taken using a simplified form of Eq. (C2) for d = 3, but 
we here generalize the calculation to arbitrary d. The 
crucial observation [13, Eq. (C21)] is that when z — > 00, 
e~ z \/2irzl n (z) -f 1. In the jamming limit the term in 
the second line of Eq. (C2) therefore disappears, be- 
cause A — > while r and u are of order 1. The re- 
maining integral can then be evaluated via the saddle 
point approximation, because both /3 and 1/A diverge. 
Consider first the case r < 1. Assuming that the sad- 
dle point u* < 1, one has to maximize the function 
— /3(1 — u) 2 — (r — u) 2 /(4A) which consistently maximizes 
u* = (4f3A + r)/(l + 40A) = (4a + rr)/(4a + r) < 1. 
Consider next the case r > 1. Assuming that in this 
case u* > 1, we have v(u) — and we thus consistently 
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Supplementary Figure 6. Growth of the isostaticity z ~ 2d (solid lines) plateau with d=3-10 for (a) HS at p — 10 10 and (b) 
SS at e ~ IO -20 . The HS-SS contact regime (c)-(d) collapses remarkably well for all d, and the G-RT predictions (red line) 
are similarly accurate as in d — 3. The plateau height also consistently decay from d=3 to 8 (insets), (e)-(f) The HS-SS 
quasi-contact power-law growth is also robustly conserved, with a constant a — 0.42(2) (blue line). The fit to Eq. (B8) is also 
provided (green line), (g) The quasi-contact coefficient C is such that the region where this regime can be observed shrinks 
with increasing d. 
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Supplementary Figure 7. (a) Fraction of rattlers in HS com- 
pressions (dashed line) and in SS energy minimizations from 
above (solid line) . The fraction of particles left outside of the 
force network, the rattlers, also vanish with increasing d. For 
both protocols, the results suggest their fraction disappears 
exponentially with d. (b) Diminishing fraction of 3-member 
force loops (triangle) and growing fraction of 4-member force 
loops (square) with d for the two protocols. The force net- 
work, which is another observable for comparing the jammed 
packings, supports their structural similarity. The length of 
the force loops is also a simple measure of structural correla- 
tions. In the mean-field high d limit these loops are expected 
to become increasingly large, as the structural correlations 
vanish. The decrease in the fraction of 3-member loops and 
the growth of the fraction of 4-member loops, accompanied 
by a constant growth of the average length of the loops, is 
consistent with this scenario. 



find u* = r > 1. Replacing these expressions for u* in 
Eq. (C2) and taking the jamming limit we obtain 



q(A,T;rr 
and Eq. (CI) reduces to 



e 4„+t 



9(l-r) 



(C3) 



S (a, t; <p) = ~ g [log(27ra) + 1] + Su q (r, (p) 



2 d -V<MGo(a,' 



(C4) 



G (a,r) 



dr r [e 4 °+ T 



which replaces [10, Eqs. (D3) and (D4)]. 

The approach to jamming from above is described by 
the small r limit [10, Appendix D.2]. In this limit we can 
consider the SS Mayer function as a small perturbation of 
the HS one and use standard liquid perturbation theory 
to write 



S liq (r,<p) = SH q s M+d2 d -V< S M / drr 



Plugging this result in Eq. (C4) we then get 



S (a,r;^) = --pog(27ra) + 1] +S 1 ? S ( ¥ >) 



+ 2 d ~\y^)d drr 
Jo 



d-l 



(C5) 



(C6) 



e 4 Q+T 
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Note en passant that the cancellation of the second term 
in Eq. (C5) with a corresponding term in Gq[q.,t) is 
not surprising, as stated in [10, Appendix D.2], but 
has a deep physical interpretation. Indeed, Eq. (C4) 
shows that the "bare" SS potential e~( 1 ~ r ) l T is mod- 
ified around jamming by the presence of m — 1 addi- 
tional replicas (with m — > 0) that "renormalizc" it to 
e -(i-r) 2 /(4a+r) ; as obtained in Eq. (C3). The crucial 
point is that the latter potential does not have a sin- 
gularity when t — y 0, ensuring a smooth crossover and 
appropriate scalings around jamming. 



2. The energy prefactor 

Starting from Eq. (C6) and repeating the calcula- 
tions of Ref. [10, Appendix D.2], we finally obtain the 
quadratic scaling of the energy as a function of Aip when 
approaching jamming from above. The general expres- 
sion for the prefactor is then easily obtained. A further 
simplification is obtained by assuming that a is small at 
the jamming point, and developing the resulting expres- 
sions in powers of \fa [13]. Doing so and optimizing over 
a and r, one finally obtains 



1 



^{p) = -d\o g ( 



d 
2 



(C7) 



Si(^) = ^[2^< S (¥>)] 2 , 
t(i P )=-Z™( V >)/(2S 1 (i P )) , 
e{ip) = [E« s (^)] 2 /(45 1 (^)) = dr 2 /(8a) . 

The second line recovers the result of Ref. [13, Eq. (77)]. 
The glass close packing point ipGCP [13], which is defined 
by the complexity Y^ s (ip) = 0, is reported in Fig. 1. The 
first line shows that yj a(</3ccp) is indeed very small, be- 
ing ~ 0.01 in d — 3 and decreasing with dimension. Lin- 
earizing the last line around </?GCP & n d using that T,^j- S (ip) 
vanishes linearly, one obtains the quadratic scaling of the 
energy and its prefactor [10]. G-RT results in Supplemen- 
tary Fig. 4 have been obtained from this procedure, using 
the Carnahan-Starling equation of state in d dimensions 
for the HS liquid [13, Eq. (82)]. 

A last remark on the energy prefactor is in order. G-RT 
results in a discrepancy between the pressure computed 
from thermodynamics and that computed from the struc- 
tural (see Refs. [13, Eq. (89)] and [10, Section VII.B]). 
This difference might have its origin in the fact that only 
two-body effective replica interactions are kept in this 
treatment. Indeed in the limit d — > oo, where this ap- 
proximation should be exact, the discrepancy disappears. 
It has also been observed in Ref. [10, Section VLB] that 
a much better agreement between theory and numerical 



data is obtained if the distance from jamming Aip is cor- 
rected to account for this discrepancy. The correction 
factor obtained from the theory corresponds to the fac- 
tor needed to impose the equality in [13, Eq. (89)], so 
rescaling Aip is equivalent to rescaling eo- The rescaled 
predictions are also reported in Supplementary Fig. 4. 



3. Scaling functions 

To complete the analysis we compute the scaling func- 
tions Z±(x) predicted by G-RT. Consider first the HS 
case, approaching jamming from below. The contact 
peak of g(r) on approaching jamming is given by [13, 
Eq. (90)] 



g{r)=g{l)A Q [2^ipg{l)^{r-l) 



— A (p— (r-1) 



(C8) 



where Aq(x) = 1 — y/irxe x [1 — erf(x)]. The validity of 
the thermodynamic relation p = 1 + 2 d ~ 1 ipg(l) is here 
assumed. As we discussed above, this relation is violated 
by the theory, but the correction is here unimportant. If 
one does not want to use this relation, it is sufficient to 
replace p — > 2 d ~ 1 pg(l), but in the end this substitution 
does not affect the prediction for the scaling function. 

Integrating Eq. (C8) using Eq. (1) we get, after chang- 
ing the variable to y — P^r{s — 1) 



Z{r) = 2dpJ^ dss' 1 - 1 A (v^is- 1 



2d- 



dy 



-1) 



d-l 



(C9) 



We notice now that p ~ 1/A</j and r — 1 ~ Aip. The 
integration is therefore over an interval of order 1. The 
first term in the integrand can be neglected because for y 
of order 1, so this term goes to 1 when Aip — > 0. Finally 
we obtain in the contact region 



Z(r) 2 



i 



e 2 " 



1 - erf 



(C10) 



where x = p(r — 1). This prediction is tested in Fig. 3 of 
the main text and in Supplementary Fig. 6. 

Next consider the SS case approaching jamming from 
above, working in the jamming limit. T — > with r = 
T/m and a = A/m. In this case the calculation starts 
from [10, Eqs. (17) and (48)]. Using bipolar coordinates 
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we can write 



g( r ) 



-pv{r 



/>oo 

/ duq(A,T;u] 
Jo 



m — 1 



e 



, ru / ru\ 
e 2A \ / ir—Id-2 — 
A — \2A 



(Cll) 



We now need to improve Eq. (C3) by considering the 
quadratic corrections around the saddle point, which for 
r < 1 leads to 



q(A,T;r) ~ e '™( 4 -+-) 



4a + rr 
r(4a + r) 



y/l + 4a/r 

(C12) 

Plugging Eq. (C12) in Eq. (Cll), dropping as before the 
last term in square brackets in Eq. (Cll), and evaluating 
the integral via the saddle point approximation including 
quadratic corrections, we find for r < 1 that 



g( r ) 



e 1; 1 + — x 



1 ( , , ( 4a 



x 9 I 1 + (r - 1) I 1 + — 



d-l 



(C13) 



Plugging this result in Eq. (1), and assuming that we 
approach jamming from above, so that r <C 4a, we obtain 



for r > 1-^ that 



z(r) = 2 &)d2<V 



4n 



1 / / 

- 1+ s-1 — 
s \ T 



(C14) 



3 -%(*-i) 2 



Changing variables to y = (1 — s)^ and using the first 
line of Eq. (C7) gives 

"1 2 

z ( r ) = d \lzr I dyil-yf^e-^ . (C15) 

(l-r)^ 

Although this result is already the desired scaling func- 
tion, we can further simplify it by noting that a is small 
to write 



Z(r) = d J— I dye' 
V 7ra J (1 _ r) 4^ 



2d 



2d 



(C16) 



where in the last line we used the last line of Eq. (C7). 
The resulting prediction 



rRT 



1 — crf(x) 



(C17) 



is tested in Fig. 3 of the main text and in Supplementary 
Fig. 6. 



